#==========================================================================
#
#  R code to draw graphs for 
#  "Careful Commitments: Democratic States and Alliance Design"
#
#       Daina Chiba (dchiba@essex.ac.uk)
#       Jesse C. Johnson (j.johnson@uky.edu)
#       Brett Ashley Leeds (leeds@rice.edu)
#
#       Last modified: 2015-04-22
#==========================================================================

rm(list=ls()) # Delete any objects in the memory

# packages ----------------------------------------------------------------

library(foreign)
library(ggplot2)
library(grid)    # for unit()
library(MASS)    # for mvrnorm()

# Load data ---------------------------------------------------------------

# k-ad data: 10 times more zeros
k10 <- read.dta("CJL_data_zero10.dta")

# Alliance observations
adata <- k10[is.na(k10 $ onlyconsul)==FALSE, ]


# S-score for alliances ---------------------------------------------------

## The code chunk below produces the histogram of interest similarity 
## shown on page 5 of the Appendix. 

g <- ggplot(data = adata, aes(min_s_2cat))
g <- g + geom_histogram()
g <- g + scale_x_continuous("Similarity of Interests",
                            limits=c(-0.5,1), breaks=c(-.5,0,.5,1))
g <- g + theme_bw()
print(g)
ggsave(file="Figure_A1.pdf",width=8,height=5)



# Figure 1: Distribution of proportion democracy --------------------------

# The code chunk below produces the histograms of our key independent 
# variable, proportion of democracies, for all the alliance observations (n=536)
# and for the defensive alliance observations (n = 260)


plot.prop.1 <- data.frame(dem.prop = adata $ dem_prop * 100, 
                          model = "All Alliances (n=536)")
plot.prop.2 <- data.frame(dem.prop = adata $ dem_prop[adata $ defense==1] * 100, 
                          model = "Defense Pacts (n=260)")
plot.prop <- rbind(plot.prop.1, plot.prop.2)

g <- ggplot(data = plot.prop, aes(dem.prop))
g <- g + facet_grid(. ~ model)
g <- g + geom_histogram(binwidth=10)
g <- g + scale_x_continuous("Percent of Democratic Members",limits=c(0,110), breaks=c(0,50,100))
g <- g + theme_bw()
print(g)
ggsave(file="Figure1.pdf",width=8,height=5)



# Simulate coefficients ---------------------------------------------------

## Load coefficient estimates from 2PMs for H1 and H2
beta.tpm1 <- as.matrix(read.table("output_2PM1beta.txt", header=T))
beta.tpm2 <- as.matrix(read.table("output_2PM2beta.txt", header=T))
## var-cov matrices
vcov.tpm1 <- as.matrix(read.table("output_2PM1vcov.txt", header=T))
vcov.tpm2 <- as.matrix(read.table("output_2PM2vcov.txt", header=T))

## Draw simulated beta from MVN(b,V)
nsims <- 1000   ## Number of resampling repetitions
set.seed(1)
simb.tpm1 <- mvrnorm(n = nsims, mu = beta.tpm1, Sigma = vcov.tpm1)
simb.tpm2 <- mvrnorm(n = nsims, mu = beta.tpm2, Sigma = vcov.tpm2)

# Set X -------------------------------------------------------------------

## Representative values of covariates

# Mean values for H1
rep1_min_s_2cat <- mean(adata $ min_s_2cat)
rep1_mem_num <- mean(adata $ mem_num)
rep1_max_threat <- mean(adata $ max_threat)
rep1_wartime <- mean(adata $ wartime)

# Mean values for H2
rep2_min_s_2cat <- mean(adata $ min_s_2cat[adata $ defense==1])
rep2_mem_num <- mean(adata $ mem_num[adata $ defense==1])
rep2_offense <- mean(adata $ offense[adata $ defense==1])
rep2_consul <- mean(adata $ consul[adata $ defense==1])

# H1
beta.tpm1
covs.t1.1o <- cbind(0,  rep1_min_s_2cat, rep1_mem_num, rep1_max_threat, rep1_wartime, 1)
covs.t1.2o <- cbind(.5, rep1_min_s_2cat, rep1_mem_num, rep1_max_threat, rep1_wartime, 1)
covs.t1.3o <- cbind(1,  rep1_min_s_2cat, rep1_mem_num, rep1_max_threat, rep1_wartime, 1)

# H2
beta.tpm2
covs.t2.1o <- cbind(0,  rep2_min_s_2cat, rep2_mem_num, rep2_offense, rep2_consul, 1)
covs.t2.2o <- cbind(.5, rep2_min_s_2cat, rep2_mem_num, rep2_offense, rep2_consul, 1)
covs.t2.3o <- cbind(1,  rep2_min_s_2cat, rep2_mem_num, rep2_offense, rep2_consul, 1)

# Calculate F(XB) ---------------------------------------------------------

# H1

prT1.out.1 <- t(plogis(covs.t1.1o %*% t(simb.tpm1[, 1:length(covs.t1.1o)])))
prT1.out.2 <- t(plogis(covs.t1.2o %*% t(simb.tpm1[, 1:length(covs.t1.2o)])))
prT1.out.3 <- t(plogis(covs.t1.3o %*% t(simb.tpm1[, 1:length(covs.t1.3o)])))

# H2
prT2.out.1 <- t(plogis(covs.t2.1o %*% t(simb.tpm2[, 1:length(covs.t2.1o)])))
prT2.out.2 <- t(plogis(covs.t2.2o %*% t(simb.tpm2[, 1:length(covs.t2.2o)])))
prT2.out.3 <- t(plogis(covs.t2.3o %*% t(simb.tpm2[, 1:length(covs.t2.3o)])))

# Figure 2 (H1) -----------------------------------------------------------

plot.mat.1 <- data.frame(prob = prT1.out.1, value = "0 %")
plot.mat.2 <- data.frame(prob = prT1.out.2, value = "50 %")
plot.mat.3 <- data.frame(prob = prT1.out.3, value = "100 %")
plot.mat <- rbind(plot.mat.1, plot.mat.2, plot.mat.3)

## manipulate the color palette
cbgFillPalette <- scale_fill_manual(values=c("white", "gray80", "black"))
y.h <- 2; y.h.2 <- 1.5; y.h.3 <- 1

mean(prT1.out.1)
mean(prT1.out.2)
mean(prT1.out.3)
ci.1 <- apply(X = prT1.out.1, MARGIN = 2, FUN = quantile, probs = c(.08, .92))
ci.2 <- apply(X = prT1.out.2, MARGIN = 2, FUN = quantile, probs = c(.08, .92))
ci.3 <- apply(X = prT1.out.3, MARGIN = 2, FUN = quantile, probs = c(.08, .92))

g <- ggplot(data = plot.mat, aes(prob, fill = value))
g <- g + geom_density(alpha = 0.8) 
g <- g + cbgFillPalette + theme_bw()
g <- g + scale_x_continuous("Predicted Probabilities of Consultation Only", limit=c(0,1),expand=c(0,0))
g <- g + scale_y_continuous("Density",limits=c(0,20))
g <- g + theme(legend.position = c(.85,.85), 
               legend.background = element_rect(fill="white",colour="black"),
               legend.title = element_text(size=12), 
               legend.text = element_text(size = 12),
               legend.key.width = unit(1,"cm"), 
               axis.ticks = element_blank(), axis.line=element_line(colour=NA), 
               axis.title.y = element_text(size=12,angle=90), axis.title.x = element_text(size=12,angle=0),
               axis.text.y = element_text(size=11), axis.text.x = element_text(size=11))
g <- g + geom_segment(aes(x = ci.1[1], y = y.h, xend = ci.1[2], yend = y.h))
g <- g + geom_segment(aes(x = ci.2[1], y = y.h.2, xend = ci.2[2], yend = y.h.2), colour="gray40")
g <- g + geom_segment(aes(x = ci.3[1], y = y.h.3, xend = ci.3[2], yend = y.h.3), colour="white")
g <- g + labs(fill = "Percent of\nDemocratic Members")
g <- g + geom_text(label = "mean = 0.20", x=0.18, y=19)
g <- g + geom_text(label = "mean = 0.35", x=0.36, y=17)
g <- g + geom_text(label = "mean = 0.55", x=0.55, y=9.5)
print(g)
ggsave(file="Figure2.pdf",width=8,height=5)


# Figure 3 (H2) -----------------------------------------------------------

plot.mat.1 <- data.frame(prob = prT2.out.1, value = "0 %")
plot.mat.2 <- data.frame(prob = prT2.out.2, value = "50 %")
plot.mat.3 <- data.frame(prob = prT2.out.3, value = "100 %")
plot.mat.d <- rbind(plot.mat.1, plot.mat.2, plot.mat.3)

ci.1 <- apply(X = prT2.out.1, MARGIN = 2, FUN = quantile, probs = c(.08, .92))
ci.2 <- apply(X = prT2.out.2, MARGIN = 2, FUN = quantile, probs = c(.08, .92))
ci.3 <- apply(X = prT2.out.3, MARGIN = 2, FUN = quantile, probs = c(.08, .92))

mean(prT2.out.1)
mean(prT2.out.2)
mean(prT2.out.3)

## ggplot: ggplot histogram 
g <- ggplot(data = plot.mat.d, aes(prob, fill = value))
g <- g + geom_density(alpha = 0.8) 
g <- g + cbgFillPalette + theme_bw()
g <- g + scale_x_continuous("Predicted Probabilities of Defense Conditionality", limit=c(0,1),expand=c(0,0))
g <- g + scale_y_continuous("Density",limits=c(0,11))
g <- g + theme(legend.position = c(.85,.85), 
               legend.background = element_rect(fill="white",colour="black"),
               legend.title = element_text(size=12), 
               legend.text = element_text(size = 12),
               legend.key.width = unit(1,"cm"), 
               axis.ticks = element_blank(), axis.line=element_line(colour=NA), 
               axis.title.y = element_text(size=12,angle=90), axis.title.x = element_text(size=12,angle=0),
               axis.text.y = element_text(size=11), axis.text.x = element_text(size=11))
g <- g + geom_segment(aes(x = ci.1[1], y = y.h, xend = ci.1[2], yend = y.h))
g <- g + geom_segment(aes(x = ci.2[1], y = y.h.2, xend = ci.2[2], yend = y.h.2), colour="gray40")
g <- g + geom_segment(aes(x = ci.3[1], y = y.h.3, xend = ci.3[2], yend = y.h.3), colour="white")
g <- g + labs(fill = "Percent of\nDemocratic Members")
g <- g + geom_text(label = "mean = 0.40", x=0.38, y=10.7)
g <- g + geom_text(label = "mean = 0.52", x=0.54, y=9.3)
g <- g + geom_text(label = "mean = 0.64", x=0.67, y=5.4)
print(g)
ggsave(file="Figure3.pdf",width=8,height=5)



## End of File